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Threshold phenomena in erosion driven by subsurface flow 

Abstract. We study channelization and slope destabilization driven by subsurface (ground- 
water) flow in a laboratory experiment. The pressure of the water entering the sandpile 
from below as well as the slope of the sandpile are varied. We present quantitative un- 
derstanding of the three modes of sediment mobilization in this experiment: surface ero- 
sion, fluidization, and slumping. The onset of erosion is controlled not only by shear stresses 
caused by surfical flows, but also hydrodynamic stresses deriving from subsurface flows. 
These additional forces require modification of the critical Shields criterion. Whereas sur- 
face flows alone can mobilize surface grains only when the water flux exceeds a thresh- 
old, subsurface flows cause this threshold to vanish at slopes steeper than a critical an- 
gle substantially smaller than the maximum angle of stability. Slopes above this criti- 
cal angle are unstable to channelization by any amount of fluid reaching the surface. 



1. Introduction 

Unlike water, a layer of sand will not flow unless its 
surface is inclined beyond a characteristic angle, known as 
the maximum angle of stability [Duran, 1999]. This simple 
fact translates into a host of threshold phenomena wherever 
granular material is found. Many such phenomena play a 
crucial role in the erosion of Earth's surface, and very likely 
manifest themselves in the richness of the patterns exhibited 
by drainage networks. 

Depending on geological, hydrological, and climatological 
properties, erosion by water is mainly driven either by over- 
land flow or subsurface flow. The former case occurs when 
the shear stress imposed by a sheet flow exceeds a thresh- 
old [Horton, 1945; Loewenherz , 1991; Dietrich et at, 1992; 
Montgomery and Dietrich, 1992; Howard, 1994]. Erosion 
in the latter case — known as seepage erosion, or sapping — 
occurs when a subsurface flow emerges on the surface. Here 
the eroding stresses derive not only from the resulting sheet 
flow but also the process of seepage itself [Iverson and Ma- 
jor, 1986]. The onset of erosion for both overland flow and 
seepage is threshold-dependent, but the additional source 
of stress in the case of seepage has the potential to create 
significantly different erosive dynamics. 

Here we study the seepage case. Whereas the case of 
Horton overland flow has been extensively studied [Smith 
and Bretherton, 1972; Dietrich et al, 1992; Montgomery 
and Dietrich, 1994; Dunne et al, 1995; Izumi and Parker, 
1995, 2000], seepage erosion has received less attention. 
Dunne [1980, 1990] suggests that erosive stresses due to 
seepage are more widespread in typical environments than 
commonly assumed. He also provides a detailed description 
of seepage erosion in the field, together with a discussion of 
the various factors that influence its occurrence. Another 
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focus of attention has been the controversial possibility that 
many erosive features on Mars appear to have resulted from 
subsurface flows [Baker, 1982; Higgms, 1982; Laity and Ma- 
lm, 1985; Baker, 1990; Aharonson et al., 2002]. Although 
the importance of seepage stresses in erosion have been re- 
alized by Howard [1988] and Howard and McLane [1988], 
comprehensive quantitative understanding is difficult to ob- 
tain. The complexity arises from the interdependent motion 
of the sediment and fluid — the "two-phase phenomenon" 
[Yalin, 1977] — which, of course, is common to all problems 
of erosion. 

To further understand seepage erosion, we proceed from 
experiments [Schumm et al., 1987]. Questions concern- 
ing the origin of ancient Martian channels have motivated 
considerable experimental work in the past [Kochel et al., 
1985; Kochel and Piper, 1986; Howard and McLane, 1988; 
Kochel et al., 1988]. The process of seepage erosion has also 
been studied as an example of drainage network develop- 
ment [Gomez and Mullen, 1992]. Our experiments, follow- 
ing those of Howard and McLane [1988]; Owoputi and Stolte 
[2001] and others, are designed to enable us to construct a 
predictive, quantitative theory. Consequently, they stress 
simplicity and completeness of information. Although our 
setup greatly simplifies much of Nature's complexity, we ex- 
pect that at least some of our conclusions will improve gen- 
eral understanding, and therefore be relevant to real, field- 
scale problems. 

A previous paper by Schdrghofer et al. [2004] provided 
a qualitative overview of the phenomenology in our experi- 
ment. It described the main modes of sediment mobilization: 
channelization, slumping, and fluidization. Here we provide 
quantitative understanding of the onset and transitions be- 
tween these modes. 

Our emphasis is on the threshold phenomena associated 
with the onset of erosion, which we will ultimately character- 
ize in the same way that others [Duran, 1999] have charac- 
terized the onset of dry granular flow beyond the maximum 
angle of stability. This involves a construction of a gen- 
eralized Shields criterion [Howard and McLane, 1988] valid 
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in the presence of seepage through an inclined surface. A 
major conclusion is that the onset of erosion driven by seep- 
age is significantly different from the onset of erosion driven 
by overland flow. We find that there is a critical slope s c , 
significantly smaller than the maximum angle of stability, 
above which the threshold disappears. Therefore any slope 
greater than s c is unstable to erosion if there is seepage 
through it. This result is similar to well-known conclusions 
for the stability to frictional failure of slopes with uniform 
seepage [Taylor, f965; Iverson and Major, 1986; Howard and 
McLane, 1988]. An important distinction in our work, how- 
ever, concerns the mode of sediment mobilization and its 
local nature. The existence of the critical slope for seepage 
erosion may provide a useful quantitative complement to the 
qualitative distinctions between seepage and overland flow 
that have already been identified [Laity and Malin, 1985]. 

The remaining modes of sediment mobilization, fluidiza- 
tion and slumping, are modeled using well established ideas 
[Iverson and Major, 1986]. The result of applying these 
ideas together with the generalized Shields criterion pro- 
vides a theoretical prediction of the outcomes of the exper- 
iment, i.e., a phase diagram. Agreement between theory 
and experiment is qualitative rather than quantitative. We 
nevertheless believe that our theoretical approach is funda- 
mentally sound and that better agreement would follow from 
improved experimental procedures. 

2. Experiment 



L 




Figure 1. Schematic of the experiment as detailed in 
Schdrghofer et al. [2004] as well as the setup for comput- 
ing the bulk and surface water flow. 

In our experimental setup, first introduced by Schdrghofer 
et al. [2004], a pile of identical cohesionless glass beads 
0.5mm in diameter is saturated with water and compacted 
to create the densest possible packing. It is then shaped 
into a trapezoidal wedge inclined at an angle 9 with slope 
s = tan# as shown in Fig. 1. The downslope length of the 
wedge is L — 90 cm, its width across the slope is 119 cm, 
and its height in the middle is approximately 11 cm. Water 
enters the sandpile underneath through a fine metal mesh 
and exits at the lower end of the pile through the same kind 
of mesh. A constant head at the inlet is maintained by 
keeping a constant water level H in the reservoir behind the 
sandbox with the help of an outflow pipe. The slope s of 
the pile and the water level H are the control parameters of 
the experiment. The degree of packing of the granular pile 
is the variable most difficult to control. 

Our particular method of feeding water into the sandpile, 
similar to that of Owoputi and Stolte [2001], can be moti- 
vated in three ways. The most important justification is the 
fact that the amount of water flowing on the surface can be 
finely controlled in our geometry. This feature is essential 
in probing the onset of erosion. Second, our setup allows us 
to access heads H larger than the height of the pile, which 



therefore allows us to explore dynamic regimes unavailable 
if water enters the pile through a mesh in the back. Third, 
a similar seepage water flow geometry can exist in the field 
wherever water travels beneath an impermeable layer that 
terminates. 

We have performed two types of experiments: steady and 
non-steady. For a fixed water level and in absence of sed- 
iment motion, water flow reaches steady state. By moni- 
toring the total water flux through the system we estimate 
the time to reach steady state to be approximately ten min- 
utes. To explore the onset of sediment motion, we raised 
the water level H in small increments, waiting each time 
for steady state to be established. Due to the particular 
shape of the bulk flow in our experiment, surface flow exists 
over a finite region of the surface. The width of this seep- 
age face and therefore the depth of the surface flow can be 
tuned by changing H . Because of the finite extent of surface 
flow, its depth and therefore the viscous shear stress reaches 
a maximum at a certain location. Thus, by increasing H 
we can continuously tune the maximum shear stress expe- 
rienced by the surface grains. The maximum shear stress 
reaches a critical value for the onset of sediment motion in 
a certain location on the slope. As we show below, we can 
compute where the maximum shear stress occurs and thus 
can reliably detect the onset of sediment motion visually 
because we focus our attention on this location. Once sed- 
iment begins to move, channels form almost immediately. 
These channels grow in length, width, and depth. An ex- 
ample of the evolving channel network is shown in Fig. 2. 
Depending on the slope, as the channels deepen, the pile 
becomes unstable to fluidization or slumping. For slopes 
lower than approximately 0.05, the fluidization threshold is 
reached before sediment is mobilized on the surface. 




Figure 2. Examples of the channel network. Top view 
of the slope of the eroding sandpile. Branched channel 
heads migrate up the slope (upward) while the eroded 
sediment is deposited in the braided wash downslope. 
The horizontal size of the image is approximately lm. 
The slope of the pile is s = 0.1 and the water level 
H = 16.15 cm. 

We also explored the non-steady evolution of the bulk and 
surface water flow and resulting sediment motion by raising 
the water level H to some higher value from zero. In this 
case one of three things can happen. The pile can be flu- 
idized within a few seconds or fail by slumping as shown in 
Fig. 3. If this does not occur, the water emerges on the sur- 
face just above the inlet. A sheet of water then washes down 
the slope of the pile. During this initial wash, sediment is 
mobilized and incipient channels form. These channels grow 
during subsequent relaxation of the bulk water flow towards 
steady state. Because of the initial wash's erosive power, 
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channels are able to form and grow for lower water pres- 



sures than in steady experiments. 
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Figure 3. Example of a slumping sandpile viewed at 
an approximately 45° angle to the slope after the water 
flow has been stopped. The width of the imaged region 
is approximately lm. Slumping happens along a convex 
upward arc which looks darker because it is deeper and 
therefore wetter. 

Outcomes of a large number of non-steady experiments 



and several steady experiments for varying slope s and the 



water level H are summarized in the phase diagram in Fig. 



4. Each symbol in the plot represents one experiment. The 



sediment is either immobile (stable seepage), or it is mobi- 



lized on the surface where channels form (channelization) 



or in the bulk (slumping or fluidization). In several ex- 



periments, slumping or fluidization happened after chan- 



nels formed and grew. In the following sections we describe 



the computations that allow us to construct the theoretical 



boundaries between the three different modes of sediment 



mobilization in our experiment. 



Figure 4. Experimental and theoretical phase diagram 
in the parameter space of slope and water level. Exper- 
iments that yielded no erosion are denoted by X; those 
that produced channels are indicated by •; and those that 
produced fluidization and/or slumping within one hour 
of the beginning of the experiment are represented by V. 
The straight line and gray-shaded curves are theoretical 
predictions for the boundaries separating the four regions 
indicated by their labels. The thickness of the lines indi- 
cates uncertainty in the theory. The boundary between 
the uneroded and channelized states is reasonably well 
approximated by our theory. The theoretical boundaries 
for fluidization and slumping, however, appear to overes- 
timate the critical water level, possibly as a result of in- 
homogeneities, dynamic changes in the sandpile's shape, 
or from the assumption of a steady state. 



3. Calculation of the water flow 

Whereas steady-state flow can be readily characterized 
quantitatively, non-steady flow characterization requires 
knowledge of the water-table dynamics. However, the the- 
ory of the water-table dynamics is less well established than 
that of the flow through the bulk of a porous medium. Also, 
our steady-state experiments probe all aspects of sediment 
dynamics. We can therefore focus on the quantitative char- 
acterization of the steady-state flow. 

To study the onset of erosion quantitatively we need to be 
able to establish a correspondence between the experimen- 
tally measurable quantities such as the slope s, the water 
level H , the size of the seepage face, and the water fluxes. 
The seepage and surface fluxes are the most difficult to mea- 
sure. In this section we set up their computation. The 
computation is designed to enable us to infer water fluxes 
indirectly by measuring the size of the seepage face. In the 
following sections we will use this computation to quantify 
the onset of erosion and to compute the slumping and liq- 
uefaction boundaries of the channelization phase diagram 
shown in Fig. 4. 

Fig. 1 specifies the key quantities and coordinate systems 
we use in computing the fluxes. The flow profile is inde- 
pendent of the ^-coordinate across the slope of the sand- 
pile except near the side walls of the box. We therefore 
treat the box as if it were infinitely wide. Flow is then two- 
dimensional and the specific discharge vector q is in the y-z 
plane. We will use two coordinate systems. As shown in Fig. 
1, the z' coordinate is measured vertically from the bottom 
of the box while the z coordinate is the normal distance 
away from the surface of the pile. The flow is governed by 
Darcy's law, 

q=-ifVV, (1) 
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where K is the scalar hydraulic conductivity, and 

1p = 7T + Z 



(2) 



ponent normal to the surface, i.e., there is either exfiltration 
or infiltration. 



is the total hydraulic head of the pore water. Both q and 
K have units of velocity while the scaled pore pressure 
7r = p/(pg) has units of length. Here p is the density of 
water and g is the magnitude of the acceleration of gravity. 
We have measured K via a U-tube relaxation experiment. 
To do so we created a water level difference A/i between the 
two arms of a transparent U-shaped tube of width I partially 
filled with glass beads. The rate of change of A/i is given by 
KAh/£. By measuring the rate of change of Ah we deduced 
the value of the hydraulic conductivity K — 3.0 ± 0.1 mm/s 
(Schdrghofer et al. [2004]). Hydraulic conductivity is sen- 
sitive to the packing of the grains and is the variable most 
difficult to control in our experiment. 

Water incompressibility implies Vq = 0, therefore yield- 
ing Laplace's equation, 



V72 

V 7T 



0. 



(3) 



To compute the pore pressure n, boundary conditions 
must be specified. The walls of the box are impenetrable. 
Therefore the discharge vector is parallel to the walls. In 
other words, the flux q± in the direction n normal to the 
walls vanishes. Thus q± = nq = 0. Because the glass beads 
in our experiment are small, capillarity is important. When 
a tube filled with glass beads is lowered into a reservoir of 
water, the porous bead-pack fully saturates in a region that 
extends above the surface of the water by a capillary rise n c . 
We measured tv c — 25 mm for our material. The capillary 
rise is a measure of the average radius of the water menisci 
at the edge of the fully saturated zone. The pore pressure at 
the edge of the fully saturated zone is — n c (without loss of 
generality we set the atmospheric pressure to zero). Water 
can rise above the fully saturated zone through the smaller 
pores and narrower throats. Thus a partially saturated cap- 
illary fringe exists above the fully saturated zone. However, 
in this fringe the water is effectively immobile since it is 
confined to the smaller pores and narrower throats. 

Since water flows only in the fully saturated zone, we 
define the water table to be at its edge. Thus, the pore 
pressure at the water table is equal to the negative capillary 
rise 7r = — 7r c . In steady state the discharge vector is par- 
allel to the water table. This extra condition allows us to 
determine the location of the water table in steady state. 

We neglect the pressure drop across the inlet mesh. 
Therefore, the pore pressure at the inlet mesh is -k = H — z . 
The boundary conditions at the surface of the sandpile and 
at the outlet mesh are more subtle. When no water seeps 
out, i.e., when the discharge vector is parallel to the sur- 
face, the curvature of the water menisci between grains can 
freely adjust so that the pressure 7r can vary between zero 
and — 7r c . Therefore when — tv c < tt < 0, no seepage occurs. 
Otherwise, the pore pressure equals the atmospheric pres- 
sure 7r = (we neglect the pressure exerted by the thin layer 
of water on the surface) , and the discharge vector has a com- 
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Figure 5. (a) Steady state location of the water table 
(dashed line) obtained by solving the Laplace's equation 
as detailed in the text. Note that the water table is al- 
most entirely at the surface of the sandpile except a for 
small region at the bottom of the sandpile. (b) Pore pres- 
sure at the water table and surface flow. Vertical scale 
for the pressure at the water table is given by the nega- 
tive capillary rise — tt c , or the pressure at the water table 
when it is below the surface of the pile. Note that seepage 
occurs only where the pore pressure reaches atmospheric 
pressure. Slope is s — 0.1, water level H = 15 cm. 



To obtain the steady state location of the water table, we 
guess its position and solve Laplace's equation (3) with the 
7r = — 7r c boundary condition on the water table. We then 
move the water table in the direction of the local discharge 
vector by an amount proportional to its length. Iteration of 
this procedure converges to the steady-state position of the 
water table. An example is shown in Figure 5. 

Once the steady flow pattern is known, we can cal- 
culate the overland water flux Q(y) by integrating the 
one-dimensional continuity condition which states that the 
downslope derivative of the overland flux is equal to the 
seepage flux: 



dy 



z q 



(4) 



y=0 



4. Onset of erosion and channel incision 

In this section we assume, based on direct observation, 
that the onset of channel incision coincides with the onset 
of erosion (i.e., we never observed a homogeneously erod- 
ing state). In other words, when the overland water flux 
becomes strong enough to carry grains, the flow of sedi- 
ment becomes immediately unstable to perturbations trans- 
verse to the downslope direction and incipient channels form 
[Smith and Bretherton, 1972]. Using this assumption and 
the calculation of the overland water flux we can deduce the 
threshold condition for the onset of erosion. 

It is universally assumed after Shields [1936] that the hy- 
drodynamic stresses exerted on the sandpile by the fluid 
flowing on its surface determine whether cohesionless gran- 
ular material is entrained. In the limit of laminar flow, the 
dominant hydrodynamic stress is the viscous shear stress Tt,. 
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Appropriately scaled this shear stress is termed the Shields 
number [Yalin, 1977], defined by 



(p s -p)gd' 



(5) 



where p B is the density of the granular material, d is the 
grain diameter (d = 0.5 mm in our experiment), and the 
surface is not inclined. 

The conventional Shields number (5) is the ratio between 
the horizontal force exerted by the flow and vertical force 
due to grain's weight. To generalize the notion of the 
Shields number to the situation with seepage through an 
inclined surface, we make two changes in Eq. (5). We first 
add the tangential component of the seepage force density 
f = —pgVip acting over a length d to the numerator of (5). 
The numerator thus becomes Tb + d(y ■ f). Note that we did 
not include the tangential component of the grain's weight to 
the numerator. Defined in this way, the generalized Shields 
number measures the effect of the fluid: both the bulk as 
well as the surface flows. 




Seepage 
Force 



Resultant 



Figure 6. The denominator of the generalized Shields 
number is the projection onto the z-axis of the resultant 
of the grain's weight and the seepage force both scaled 
by yrd 2 /6. 



Second, we replace the denominator of Eq. (5) by the 
resultant (vectorial sum) of the seepage force on a grain and 
its submerged weight, as shown in Fig. 6, both scaled by 
7rd 2 /6 (d 2 to obtain stress as in the numerator and n/6 for 
agreement with the conventional Shields number), projected 
onto the z-axis. According to Martin and Aral [1971] the 
grains on the surface of the bed experience a seepage force 
roughly half as large as the grains several layers deep. Con- 
sequently, we assume that the seepage force is reduced by a 
factor of a w 0.5; therefore 



r b - apgd (dip/dy) 



(p s — p)gd cos 9 + apgd (dip/dz)' 



(6) 



where 9 is the inclination angle of the surface. The impor- 
tance of the seepage stresses for the criterion for the onset 
of erosion was previously realized by Howard and McLane 
[1988]. It can be shown that their equation (10) expressing 
marginal stability of a surface grain is equivalent to writing 
r* = tancp, the tangent of the angle of internal friction. 

The generalized Shields number Eq. (6) is a measure of 
the relative importance of the tangential and normal forces 
acting on a grain at the surface of the sandpile. There- 
fore, we expect r* to be a control parameter for erosion. In 
other words, there exists a critical Shields number r * , such 
that when r* < r*, surface grains are immobile, and when 
r* > r*, sediment is mobilized. Note that (6) reduces to the 
classical definition of the Shields number for a flat surface 
without seepage. Also note that since we did not include the 



tangential component of grain's weight, the critical Shields 
number at the onset of sediment motion vanishes when the 
inclination angle 9 reaches the maximum angle of stability. 

Although we obtain the seepage force density f as a re- 
sult of computing the pore- water pressure, to calculate the 
boundary shear stress r b we must estimate the thickness of 
the surface water layer h{y). Since this thickness changes 
slowly in the downslope direction, we can approximate the 
surface flow by the steady flow of a uniform layer of viscous 
fluid. Also, the surface water flux is small enough for tur- 
bulence to be of no importance. The thickness h of laminar 
surface flow for a given flux Q is [Landau and Lifshitz, 1987] 



h 



( 3Q77 A 1/3 
V P9 sin t 



(7) 



where 7/ is the viscosity, while the viscous shear stress ex- 
erted on the sandpile is 



Tb = pgh sin 9. 



(8) 



The particle Reynolds number can then be calculated using 
the bottom shear stress (8) and shear velocity v* — y/Th/p 



Re = 



v*d 



(9) 



where v = r//p is the kinematic viscosity of water. We esti- 
mate that in our experiments, this particle Reynolds num- 
ber varies between 5 and 20 depending on the slope s of the 
pile and the water level H . We verify this estimate of the 
Reynolds number by a direct measurement of the thickness 
of the surface flow. We find that this thickness is several 
grain diameters. This justifies the laminar flow assumption 
used in obtaining Eq. (7). 

Using (8), the Shields number (6) can now be conve- 
niently rewritten as 



(h/d) sin 9- a (dip/dy) 
(p s /p — 1) cos 9 + a (dip/dz) ' 



(10) 



This expression can be further simplified by noting that 
-K — along the seepage face. Therefore dtp/dy = —sin 9 
at the surface wherever there is overland flow. We arrive at 
the final expression for the modified Shields number which 
depends on the surface flow thickness h, the normal compo- 
nent dip/dz of the seepage force density at the surface, and 
the seepage force reduction factor a 



(h/d + a) sin 9 



(p a /p- l)cos9 + a (dip/dz)' 



(11) 



In our geometry, both the surface and the seepage water 
fluxes reach a maximum somewhere along the slope. There- 
fore the Shields number has a maximum value as well. Below 
we calculate this maximum Shields number in steady state 
for a given slope s and water level H. 

5. Critical slope for the onset of seepage 
erosion 
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Figure 8. Schematic of the pore pressure n as a func- 
tion of the downslope coordinate y at the onset of seepage 
H — H B and just above. 



Figure 7. Computed maximum Shields number ver- 
sus the water level H for a pile of slope s — 0.1. At 
H — H B water first seeps through the surface and the 
Shields number jumps to a nonzero value. Afterwards 
it increases rapidly as (H — Hs) 1 / 3 (solid line). Inset: 
corresponding size of the computed seepage face. 

We now explore the consequences of seepage for the phe- 
nomenology of the onset of erosion. Because of the addi- 
tional force on the surface grains, seepage flow is more ero- 
sive than overland flow. This notion is reflected quantita- 
tively in the generalized Shields number. Let us examine 
how the maximum Shields number varies with the water 
level H in our experiment. A representative graph of the 
maximum Shields number versus the water level is shown in 
Fig. 7. Below a water level H B (s) that is a function of the 
slope, no water seeps out to the surface of the pile. Even 
though the water table may be at the surface, the pressure 
at the water table is below atmospheric pressure and cap- 
illarity prevents seepage. When H — H B , i.e., exactly at 
the onset of seepage, the pressure reaches it = at some 
point on the surface. Since the seepage flux is still zero, 
dtp/dz = along the wet part of the surface. Therefore, 
just above the seepage onset, when the water layer thick- 
ness and the seepage flux are both infinitesimally small, the 
maximum Shields number is 

™-1&=t < I2 > 

In contrast to overland flow, the consequence of seepage is 
that as soon as the water emerges on the surface, the maxi- 
mum Shields number is some non-zero value which depends 
on the slope. This also implies that there exists a critical 
slope s c such that r s (s c ) is equal to the critical Shields num- 
ber t* , i.e., 



s c = r c *(p s /p - I) /a. (13) 



For slopes greater than s c seepage is always erosive. Note 
that for low-density particles this critical slope can be arbi- 
trarily small. The expression for the critical slope for seep- 
age erosion in Eq. (13) is analogous to well-known formulas 
for stability of slopes to Coulomb failure due to uniform 
seepage [Iverson and Major, 1986]. Our result applies lo- 
cally to the point where non-uniform seepage first emerges 
on the surface. In this situation, the pile is generally stable 



We now show that above H B , the maximum Shields num- 
ber increases rapidly as a 1/3 power of the water level excess 
(H — ff s ) 1//3 . At the onset of seepage, i.e., when H = H B , 
the pressure at the water table reaches atmospheric pres- 
sure -K — at some point P located at z — and y — y B on 
the surface. Even though the water table is at the surface, 
there is no seepage anywhere, i.e., dtp/dz — 0. Because the 
pressure is smooth, it can be approximated by a quadratic 
function near this point so that tv « bz 2 — c(y — y B ) 2 , where 
b and c are constants with appropriate dimensions. When 
the water level is raised by a small increment SH = H — H B , 
the lowest order change in the head at the water table is an 
increase of SH with the exception of the region where this 
increase would lead to a positive pressure. As illustrated in 
Fig. 8, in this region the pore pressure n is set to 0, and thus 
this region becomes the seepage face. The width of the seep- 
age face W scales like the square root of SH, i.e., W ~ SH 1 / 2 
as seen in the inset of Fig. 7. The seepage flux can be es- 
timated by noting that the hydraulic head is modified by 
an amount SH over a vertical region of order W . Therefore 
we obtain dtp/dz ~ SH/W ~ SH 1/2 . The total surface flux 
therefore scales like the product of the seepage flux and the 
width W of the seepage face, i.e., Q ~ W dtp/dz ~ SH. The 
lowest order change in the maximum Shields number is due 
to the change of the surface flow depth h ~ Q 1 ^ 3 ~ SH 1 / 3 . 
Thus as we claimed above, just above the water level H B for 
the onset of seepage, 

r* «t,*(*) + a( ( )(ff-fl,) 1/s , (14) 

where the constant a is a function of the slope. Variation 
with water level of the computed maximum Shields number 
shown in Fig. 7 is consistent with Eq. (14). 

6. Measurement of the critical Shields 
number 

In the previous sections we have detailed the way of cal- 
culating the bulk and surface water fluxes in our experiment 
and the resulting maximum generalized Shields number. In 
this section we use this calculation to examine the onset of 
the sediment flow and channelization. Our first goal is to 
measure the threshold or critical Shields number required 
for the mobilization of sediment. We then use this measured 
value of the critical Shields number to predict the outcome 
of steady-state experiments for various values of the slope 
and the water level and thus compute the channelization 
boundary in the phase diagram in Fig. 4. 

The actual maximum Shields number in the experiment 
differs from the quantity calculated in Eq. (11). In addition 
to random errors in the measurements of the pile dimensions 
and water level, there are several sources of systematic error. 
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For example, the pressure drop across the inlet mesh results 
in a lower effective hydraulic head. Also, our measurement 
of the capillary rise n c is dependent on a visual estimate of 
the fully saturated zone and thus can be a source of sys- 
tematic error. We indeed find that the size of the seepage 
face calculated for a particular water level H is greater than 
measured in the experiment. However, the size of the seep- 
age face translates directly into the surface water flux and 
therefore the maximum Shields number. 




0.02 1 1 1 1 1 

20 40 60 80 

Size of seepage face, cm 

Figure 9. Computed maximum Shields number versus 
the seepage face size for three different slopes. We use 
this calculation to infer the Shields number from the ex- 
perimentally measured size of the seepage face. 

The inset of Figure 7 shows the typical dependence of 
the size of the seepage face on the water level. The vari- 
ation of the maximum Shields number with the size of the 
seepage face is shown in Fig. 9 for three different slopes. We 
use this computed correspondence between the size of the 
seepage face and the maximum Shields number to infer the 
maximum Shields number in the experiment by measuring 
the size of the seepage face. 

To measure the critical Shields number we raise the water 
level H in increments of a few millimeters at a time. Each 
time the water level is increased, the seepage flow is allowed 
to reach a steady state. In each of these steady states we 
measure the seepage face size and infer the corresponding 
maximum Shields number. Eventually, sediment is mobi- 
lized and we record the size of the seepage face and compute 
the corresponding maximum Shields number. This number 
is an upper bound on the critical Shields number for our 
granular material at that particular slope. The lower bound 
on the critical Shields number is obtained from the largest 
seepage face at which no sediment is moving or sediment mo- 
tion is only transient. Averaging over several experiments 
with slope s = 0.1 we estimate the critical Shields number 
to be 

t* = 0.085 ± 0.005. (15) 

It is not obvious that the generalized critical Shields num- 
ber for the onset of seepage driven erosion should coincide 
with the critical Shields number for overland flow. How- 
ever our measured value of the critical generalized Shields 
number is within the scatter of the existing data for over- 
land flow summarized in Buffington and Montgomery [1997]. 
Our measurement the critical generalized Shields number is 
equivalent to measuring the angle of internal friction due to 
the correspondence of our definition of r* (6) and Howard 
and McLane's equation (10). 

Deviations from flatness of the pile's surface result in the 
fluctuations of the thickness of the surface water film. As a 



result, the maximum bottom shear stress in the experiment 
is systematically greater than that calculated at a given size 
of the seepage face. Thus the Shields number calculated for 
a particular size of the seepage face is the lower bound on 
the actual Shields number in the experiment. 

In principle, the critical Shields number should vary with 
the slope of the pile. Evidence for this is the fact that at 
the maximum angle of stability any additional forcing from 
the water flowing over the bed mobilizes sediment. Since 
it is reasonable to assume that the critical Shields number 
is continuous and monotonic, we arrive at the notion that 
it decreases monotonically with slope and vanishes at the 
maximum angle of stability. For small slopes the critical 
Shields number is expected to decrease as the cosine of the 
inclination angle since this is the lowest order change in the 
stabilizing effect of gravity. For most slopes in our exper- 
iments, cos 8 is within a few percent of unity and thus we 
can ignore the variation of the critical Shields number with 
slope. 

This assumption allows us to predict the water level 
at which erosion and therefore channelization should com- 
mence in our experiment. In Fig. 4 a boundary is drawn 
between regions where sediment is expected to mobilize and 
remain immobile. To obtain this line we computed for each 
slope the water level H c (s) at which the Shields number is 
equal to the critical Shields number (15). Below this water 
level, i.e., when H < H c (s), the maximum Shields number 
is below critical and thus sediment is immobile. Conversely, 
for H > H c (s), the maximum Shields number is above crit- 
ical and thus sediment is mobilized and channels form. The 
channelization boundary is widened because of the uncer- 
tainty in the critical Shields number. 

Qualitative agreement of the channelization boundary 
with experiments is perhaps due to the opposite action of 
two effects. First, channelization occurs for lower water lev- 
els in non-steady experiments. This happens because in 
non-steady experiments the maximum Shields number over- 
shoots its steady-state value. The overshoot is greatest for 
small slopes. Second, a pressure drop across the inlet mesh 
and the compacted region of sand close to it has an oppo- 
site effect which increases the water level needed for chan- 
nelization. These two effects, though small, could together 
affect the accuracy of our predictions. Since these effects 
act in opposite ways, our the predictions of the calculated 
channelization water level H c (s) agree qualitatively with the 
experiments. 

7. Fluidization and slumping 

Having computed the channelization boundary in the 
phase diagram, we now pursue a quantitative description 
of the other two modes of sediment mobilization exhibited 
by our sandpile. Higher water pressures can cause the sand- 
pile to fail in one of two ways. First, an upward seepage 
force can lift sand and result in a fluidization or quicksand 
instability. Second, the pile can become unstable to slip- 
ping, slumping, or sliding. Both failure mechanisms have 
been discussed by a number of studies, e.g., those of Iverson 
and Major [1986] or Martin and Aral [1971]. 
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Figure 10. To predict fluidization, we calculate the to- 
tal seepage force on a slice of width A and height I and 
compare it with its weight. 

Fluidization occurs when at some point P in the sandpile 
the pore pressure is larger than the total hydrostatic pres- 
sure due to the weight of the sand and the water above P. 
To see this we compute the total seepage force acting on a 
slice of sand of width A between point P and point O on the 
surface of the pile directly above P. The vertical component 
of this force is (see Fig. 10) 

F v = pgA(4> P - ipo) = pc/A(7T P - no - I), (16) 

where I is the height of the slice. When this force exceeds 
the submerged weight A£(p t — p)g of the slice, the slice is 
lifted and the bed is fluidized. Here p t is the total density 
of the saturated sand, which for our sand is approximately 
2 g/cm 3 . Thus fluidization occurs when there exists points 
P and O on the surface directly above P such that 

vrp - 7T > Zpt/p. (17) 



Second, we compute the instability of an uneroded pile, 
whereas in most experiments, the pile failed after erosion 
had changed the shape of the pile. The decrease of pile's 
height due to erosion increases the head gradient in the bulk 
and thus makes the pile more prone to slumping and/or flu- 
idization. 
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Figure 11. Effective inclination angle at the surface due 
to seepage forces for s — 0.18 and three different water 
levels. For the highest water level, a region on the surface 
has an effective angle above the maximum angle of stabil- 
ity and thus the slope is unstable to slumping. The inset 
shows the plateau value of the effective inclination angle 
as a function of slope. When the plateau value reaches 
the maximum angle of stability, even a small amount of 
seepage destabilizes the pile to slumping. 



For uniform seepage this condition is equivalent to those in 
Iverson and Major [1986] and Martin and Aral [1971]. 

To construct the fluidization boundary in the phase di- 
agram (Fig. 4), we find the water level //fluid above which 
there exists at least one point in the pile for which condi- 
tion (17) is satisfied. Below this fluidization water level this 
condition is not satisfied for any point in the pile. 

In addition to fluidization the sandpile can fail by slump- 
ing. This can happen in one of two ways. Frictional failure 
can occur in the bulk of the pile due to the seepage stresses. 
Alternatively, surface avalanching can occur. To establish 
an upper bound on the water level at which the sandpile 
slumps via either mechanism we use the criterion developed 
by Iverson and Major [1986] for determining when a slope is 
destabilized by uniform groundwater seepage. Essentially it 
requires calculating the vectorial sum of the seepage and 
gravity forces acting on a small element of soil near the 
surface. When the angle between this total force and the 
downward normal to the surface, which we will call the ef- 
fective inclination angle, exceeds the maximum angle of sta- 
bility, the surface grains are destabilized. We measured the 
maximum angle of stability to be ip « 23° for dry glass 
beads. The slumping boundary in the phase diagram (Fig. 
4) is constructed by computing the effective inclination an- 
gle along the surface of the pile and noting the water level 
//slump, at which the effective inclination angle reaches the 
maximum angle of stability ip at some point of the surface. 

Figure 4 shows the critical water level at which fluidiza- 
tion and slumping should occur according to the criteria 
above. Failure occurs at systematically lower water levels in 
the experiment. There are several effects which can account 
for this difference. First, any irregularities in the construc- 
tion of the pile such as voids or surface height fluctuations 
make the pile more unstable to fluidization and slumping. 



At s « 0.2, a jump in the slumping water level is observed 
in both the experiment and the model. This jump is a purely 
geometric effect. Slumping occurs when, somewhere along 
the slope, the effective inclination angle, which includes the 
effect of the seepage force, exceeds the maximum angle of 
stability. As shown in Fig. 11, for slopes smaller than 0.2, 
the effective inclination angle is flat and develops a peak un- 
der the water inlet as the water pressure is increased. When 
the top of this peak crosses the value of the maximum angle 
of stability, the pile slumps. When the slope exceeds 0.2, 
however, the value of the plateau in the effective inclination 
angle is above the maximum angle of stability. Therefore, 
for these slopes, the pile will be unstable to slumping as soon 
as the water emerges on the surface. 

8. Conclusions 

This article reports on our progress in understanding 
seepage erosion of a simple non-cohesive granular material 
in a laboratory-scale experiment introduced in Schdrghofer 
et al. [2004]. Our ultimate goal is to construct a quantita- 
tive predictive theory of the onset and growth of the channel 
network observed in this experiment. This goal requires a 
complete sediment transport model as well as the calcula- 
tion of the relevant water fluxes. Here we obtain the latter 
and focus on the onset of erosion. 

Prediction of the onset of erosion based on the gener- 
alized Shields conjecture explains qualitatively the chan- 
nelization boundary in the experimental phase diagram. 
By invoking well established simple ideas we also roughly 
explain the fluidization and slumping boundaries in the 
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phase diagram. Greater discrepancy with the experiment 
for these boundaries indicates that better understanding of 
the slumping/fluidization mechanisms particular to our ex- 
periment is needed. 

The central result of our exploration is the introduction 
of the generalized Shields criterion for seepage erosion. As 
a consequence of seepage forces on the surface grains, the 
threshold for the onset of erosion driven by seepage is slope 
dependent. The threshold disappears at a critical slope 
s c determined by the critical Shields number for overland 
flow and the density contrast between the granular material 
and water. In most cases this critical slope is significantly 
smaller than the maximum angle of stability. We find, there- 
fore, that slopes above this critical slope are unstable to any 
amount of seepage. As a consequence, slopes that sustain 
seepage must be inclined at an angle smaller than the critical 
angle for seepage erosion. This behavior contrasts strongly 
with the threshold phenomena in erosion by overland flow, 
and therefore provides a mechanistic foundation for distin- 
guishing the two types of erosion. 
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